[ <--- prev -- ]  [ HOME ]  [ -- next ---> ]

[ full index ]


16 Special source: cosmic rays

Cosmic ray calculations can be done with FLUKA using the input commands GCR-SPE (for initialisation purposes) and SPECSOUR. In addition, several auxiliary stand-alone programs need to be used to prepare the geometry and material cards to be inserted into the input file.

Command SPECSOUR does not define only cosmic ray sources, but also a number of other pre-defined complex sources that cannot be described by the simple keywords BEAM, BEAMPOSit, BEAMAXES and HI-PROPErt. To avoid confusion, SPECSOUR input related to cosmic rays is described separately in this Chapter, in 16.7}.

A complete calculation to determine particle fluxes in the atmosphere requires:

The following options are available concerning the simulation of cosmic ray interactions in FLUKA:

The DMPJET and the superposition model can also be used together, by setting the respective energy ranges with the PHYSICS card.

A number of tools and packages have been developed for the FLUKA environment to simulate the production of secondary particles by primary cosmic rays interacting with the Earth's atmosphere. These tools, in different stand-alone versions, have already been successfully used for fundamental physics research [Bat00,Bat03a]. '

The set of FLUKA tools for cosmic ray simulation includes a set of core routines to manage event generation, geomagnetic effects and particle scoring, and the following stand-alone data files and programs :

16.1 Primary spectrum

The Galactic Cosmic Ray (GCR) component of the cosmic ray flux can be simulated up to 30 TeV/nucleon (or 500 TeV/n when DPMJET is linked). Two options are available: with the first one the actual ion composition of the flux is used (All-Particle Spectrum), while with the second option the primary flux is treated as a sum of nucleons (All-Nucleon Spectrum).


16.1.1 The All-Particle Spectrum

The ion composition of the galactic flux is derived from a code [Bad96] which considers all elemental groups from Z = 1 to Z = 28. The spectrum is modified to follow recent data sets (AMS [Alc00,Alc00a] and BESS [San00] data of 1998) up to 100 GeV according to the so-called ICRC2001 fit [Gai01]. The spectrum components are written into 28 files. The name of the files has the form (Z+phi+<PhiMV>+.spc). The first two characters of each file name are the atomic number of a different primary spectrum ion (e.g. 01:protons, 02:alpha...). They are followed by the solar modulation parameter used for generating the spectrum (7 characters) and by an extension ".spc". The ".spc" files are spectra without geomagnetic cutoff. The ".spc" files are used together with an analytical calculation of the rigidity cutoff, according to a centered dipole approximation of the Earth geomagnetic field, adapted to result in the vertical cutoff inserted into the input file (SPECSOUR command, SDUM=GCR-IONF, WHAT(2) of the continuation card), at the geomagnetic latitude and longitude of interest.


16.1.2 The High Energy All-Nucleon Spectrum

The All-Nucleon Spectrum is obtained modifying the fit of the All-Nucleon flux proposed by the Bartol group [Agr96], using the All-Particle Spectrum (16.1.1}) up to 100 GeV and data published in ICRC 2003. Fluxes are read from a file named ällnucok.dat" giving the total energy (GeV), the fluxes (E.dN/dE) and the neutron/proton ratios. This option ("All Nucleon Flux") is chosen with command SPECSOUR and SDUM = GCR-ALLF (see details in 16.7}). The user can decide whether to sample neutrons and protons from the file and to transport them using the superposition model, or to consider all neutrons as being bound in alpha particles and to transport protons and alphas. This latter choice has the advantage of taking better into account the magnetic field, which has no effect on the neutrons.

For the proton component at energies larger than 100 GeV, using the normalization obtained at 100 GeV, a spectral index gamma = -2.71 is assumed. A spectral index gamma = -3.11 is assumed above the knee at 3000 TeV. For what concerns the He component, gamma = -2.59 is used above 100 GeV and a charge-dependent knee is assumed according to the rule: E_nucleon = Z * 3000 TeV/A. Higher Z components have been grouped in CNO, MgSi and Fe sets and treated using an All-Particle spectrum with the above mentioned charge-dependent knee parameterisation.

16.2 Solar modulation

The deviation from the power law, observed below 10 GeV, is a consequence of the influence of the solar wind called solar modulation [Gle68]. Flux intensity in this energy range is anti-correlated to the solar activity and follows the sun-spot 11-year cycle. The correlation between the solar activity and the modulation of the cosmic rays flux has been studied by monitoring the flux of atmospheric neutrons. In fact, a flux of low energy neutrons (E ~ 1.E8-1.E9 eV) is produced in the interaction of primary CRs with the atmosphere and it is mostly due to low energy primaries (1-20 GeV), due to the rapid fall of the primary flux intensity with energy. One assumes that far from the solar system there exists an unmodified flux called Local Interstellar Spectrum, which is modified within the solar system by the interaction with the solar wind. This interaction is well described by the Fokker-Planck diffusion equation. Describing the solar wind by a set of magnetic irregularities, and considering these irregularities as perfect elastic scattering centres, one obtains the Fokker-Planck diffusion equation. For energies above 100 MeV this equation can be solved using the "Force Field Approximation" [Cab04]. According to this approximation, at a given distance from the Sun, for example at 1 a.u., the population of CRs at energy E_interstellar is shifted at the energy E_0 as in an energy loss mechanism due to a potential V:

 E_0 = E_interstellar + Z . V_solarwind(t)

The solar wind potential at a given distance from the Sun depends on only one parameter, the time: V = V(t). So it doesn't matter what the interstellar flux is: given a flux on the Earth at a time t, one can find the flux at another time just from the relative variation of the solar wind potential Phi. This variation can be derived from the neutron monitor counts [Bad96]. In the case of the fit used by FLUKA, an offline code [Bad96] makes use of an algorithm which takes into account a specific Phi value, or the counting rate of the CLIMAX neutron monitor [CLIMAX] to provide the prediction for the flux at a specific date or for a given value of the potential which expresses the effect of the interplanetary modulation of the local interstellar spectrum. Even if the model is not a description of the processes and of the manner in which they occur, it reasonably predicts the GCR modulation at Earth.

16.3 Atmospheric model: geometry


16.3.1 Earth atmosphere model

The FLUKA package makes use of a density vs. height profile of atmosphere. An external program containing a functional fit to this profile has been used to generate at the same time an input geometry file, together with the data cards for material description (each atmospheric layer, having its proper density, needs to be assigned a different FLUKA material). The geometry produced, and distributed with the name atmogeo.cards is a spherical representation of the whole Earth atmosphere. The material definitions and assigment contained in the file atmomat.cards correspond to the density profile of the U.S. Standard atmosphere. The cards contained in atmomat.cards shall be included by the user in her/his input file. In addition, the user can specialize this geometry to a given geomagnetic latitude and longitude with the help of the atmloc_2011.f auxiliary program. In this way, the geometry will contain only a slice of the atmosphere, centered on the given position. The local geometry file produced by atmloc_2011.f} is named atmloc.geo. The user shall rename this geometry file for further use. More auxiliary files are produced by atmloc_2011.f: the file atmlocmat.cards contain additional material assignments to be included in the input together with the ones from atmogeo.cards; the file atmloc.sur contains data used by FLUKA runtime, and normalization areas.


16.3.2 Local atmosphere model

The geometry is built using two truncated cones (TRC) whose vertex is in the centre of the Earth, the base is out of the atmosphere and the altitude (considering a geographical location in the northern hemisphere) is in the direction of the Earth radius which passes through the North Pole. The angular span between the two cones contains the atmosphere of interest for the latitude of interest. In addition there is a third cone placed in the opposite direction: its vertex is where the other two cones have the base, its base is out of the atmosphere and its height is in the direction of the Earth radius which passes through the South Pole. A similar geometry can be built for a requested latitude in the southern emisphere. So the complete geometry of the local model, built with the auxiliary program atmloc_2011.f, is made of:

16.4 Atmospheric model: density

The atmosphere can be roughly characterized as the region from sea level to about 1000 km altitude around the globe, where neutral gases can be detected. Below 50 km the atmosphere can be assumed to be homogeneously mixed and can be treated as a perfect gas. Above 80 km the hydrostatic equilibrium gradually breaks down as diffusion and vertical transport become important.

The following Table shows the U.S. Standard Atmosphere depth vs altitude and vs FLUKA atmospheric layer.